close all;
clear;
clc;

%% 1. 读取图像并二值化
img = imread("Fig01.tif");
img = img(:,:,1);
figure("Position",[185 60 1140 453]);
subplot(1,2,1);
imshow(img);
subplot(1,2,2);
histogram(img);

% 自适应二值化
binary_img = imbinarize(img, "adaptive", "ForegroundPolarity", "bright", "Sensitivity", 0.4);
figure("Position",[185 60 1140 453]);
subplot(1,4,1);
imshow(binary_img);
title("自适应二值化图像");

%% 2. 开运算（去噪声）

% 计算平均短轴长度
stats = regionprops(binary_img, 'MajorAxisLength', 'MinorAxisLength');
avgMinorAxis = mean([stats.MinorAxisLength]);
radius_open = floor(avgMinorAxis / 4); % 开运算半径
binary_img = imopen(binary_img, strel("disk", radius_open));
subplot(1,4,2);
imshow(binary_img);
title("开运算去噪声结果");
%% 3. 闭运算（连接粘连）
radius_close = ceil(avgMinorAxis / 3); % 闭运算半径
binary_img = imclose(binary_img, strel("disk", radius_close));
subplot(1,4,3);
imshow(binary_img);
title("闭运算平滑轮廓结果");
%% 4. 分水岭分割处理粘连区域
% 距离变换图
distance = -bwdist(~binary_img);
distance = imgaussfilt(distance, 0.5); % 平滑距离图
distance(~binary_img) = -Inf;

% 提取种子点
markers = imextendedmin(distance, 15);
%markers = imdilate(markers, strel('disk', 10)); % 种子点膨胀

% 限制分水岭敏感性
distance_filtered = imhmin(distance, 0.44);

% 分水岭分割
L = watershed(distance_filtered);
binary_img(L == 0) = 0; % 分割线作为背景

% 清理过小区域
cc = bwconncomp(binary_img);
stats = regionprops(cc, 'Area');
min_area = 20; % 设置最小米粒面积
for i = 1:cc.NumObjects
    if stats(i).Area < min_area
        binary_img(cc.PixelIdxList{i}) = 0; % 移除过小区域
    end
end

% 可视化结果
subplot(1,4,4);
imshow(binary_img, [0, 1]);
title('优化后的分水岭分割结果');

% 统计目标数量和面积
cc = bwconncomp(binary_img);
stats = regionprops(cc);
areas = [stats.Area];
disp(['米粒的总数量是：', num2str(length(areas))]);
disp('每个米粒的面积如下（像素）：');
disp(areas);

% 显示标记结果
label = labelmatrix(cc);
RGB_label = label2rgb(label, 'jet', 'k', 'shuffle');
figure("Position",[185 60 1140 453]);
subplot(1,2,1);
imshow(RGB_label);
title('标记的米粒区域');

% 面积直方图
subplot(1,2,2);
histogram(areas);
title('米粒面积分布');
xlabel('面积（像素）');
ylabel('频数');
